import math as m
import os


import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import h5py
import ast
from sympy import *

import copy

from sympy import *

# np.set_printoptions(suppress=True)
import copy

def draw_Logistic_bifurcation_last():
    '''
    Logistic_bifurcation需要和draw_logistic_bifurcation联用
    :return:
    '''
    init_value = 0.32
    _xlist = np.zeros(100)
    _xlist[0] = init_value
    for u_index in np.arange(0, 4.1, 0.1):
        for index in np.arange(0, 100 - 1):
            _xlist[index + 1] = u_index * _xlist[index] * (1 - _xlist[index])
            filename = r'logistic_' + str(u_index)
            with open(r'../experiment/logistic_diff/' + filename, 'w+') as file:
                for line in _xlist:
                    file.write(str(line) + '\n')
    draw_Logistic_biff()


def draw_Logistic_biff(data_pth=r'..\experiment\Logistics.txt'):
    '''
    :param data_pth: Logistic的数据路径.
    :return: null
    '''
    # 1. 装载数据
    numbers = np.loadtxt(data_pth)
    # 2. 数字转换
    numbers = np.float64(numbers)
    # 3. 画图
    plt.plot(numbers, 'r.')

    # 4. 设置坐标系名称
    plt.xlabel('times')
    plt.ylabel('values')
    plt.show()


def draw_logistic_bifurcation(data_path=r'../experiment/logistic_diff/'):
    '''
    画Logistic的分叉图
    :param data_path: 数据路径
    :return:  null
    '''
    # 设置参数
    u_min = 0
    u_max = 4.1
    u_step = 0.1
    u_value = np.arange(u_min, u_max, u_step)
    u_list = []
    value_list = np.zeros(1)
    times = 100
    # 读取迭代值
    for filename in os.listdir(data_path):
        # 1. 装载数据
        numbers = np.loadtxt(data_path + filename)
        # 2. 读取所有文件到一个变量
        numbers = np.float64(numbers)
        value_list = np.append(value_list, numbers)

    # 将控制参数扩展，构建二维坐标系
    for u_index in range(0, u_value.size):
        u = [u_value[u_index]] * times
        u_list.extend(u)

    # np化
    u_list = np.array(u_list)

    # 画图
    plt.plot(u_list, value_list[1:value_list.size], 'r.')
    plt.xlabel('u')
    plt.ylabel('value')
    plt.show()


# ##################另一种方式实现###############################
#     xlist = []
#     ulist = []
#     x0 = 0.5
#     times = 1000
#
#     for u in range(0, 401):
#         u = u / 100
#         x = u * x0 * (1 - x0)
#         for i in range(times):
#             x = u * x * (1 - x)
#             xlist.append(x)
#
#         u = [u] * times
#         ulist.extend(u)
#
#     x = np.array(ulist)
#     y = np.array(xlist)
#
#     # plt.plot(x, y, c='yellow', label='Logistic')
#     plt.scatter(x, y, s=1, alpha=0.4, label='Logistic')
#     plt.ylim((0, 1.05))
#     plt.xlim((0, 4.05))
#     plt.show()
def draw_all_bifurcation(parameters=4, times=100, step=0.1, data_path=r'../experiment/logistic_diff/'):
    '''
    可以画所有1d混沌系统参数,数据结构[参数,混沌值]，只能画2d的图，三维的三点这个函数画不了
    :param parameters:  系统参数,也是坐标系的横坐标
    :param times:  循环次数,迭代次数
    :param data_path:  数据路径
    :param step: 默认步长0.1
    :return: null
    '''
    # 设置控制参数
    parameters_min = 0
    parameters_max = parameters + 0.1
    u_step = 0.1
    u_value = np.arange(parameters_min, parameters_max, u_step)
    u_list = []
    value_list = np.zeros(1)

    # 读取所有的迭代值
    for filename in os.listdir(data_path):
        # 1. 装载数据
        numbers = np.loadtxt(data_path + filename)
        # 2. 读取所有文件到一个变量
        numbers = np.float64(numbers)
        value_list = np.append(value_list, numbers)

    # 生成控制参数控制范围，因为python里需要所有的坐标系长度是一样的
    for u_index in range(0, u_value.size):
        u = [u_value[u_index]] * times
        u_list.extend(u)
    # np化
    u_list = np.array(u_list)

    # 画图
    plt.plot(u_list, value_list[1:value_list.size], 'r.')
    # 设置坐标系
    plt.xlabel('u')
    plt.ylabel('value')
    # 显示图像
    plt.show()


def draw_bifurcate_2d_sine_logistic(times, x_init, y_init):
    '''
    draw_bifurcate_2d_sine_logistic 是画sine_logistic混沌映射的三维散点图，
    x，y坐标是两个控制参数 u，w。z坐标是迭代值，迭代值可以是x，y两种迭代值，这个x，y是迭代结果不是坐标系。
    :param times: 迭代的次数
    :param x_init: x的初始值
    :param y_init: y的初始值
    :return: 返回值为null
    '''
    # 判断当前文件夹下是否存在结果，如果有的话，需要删除。
    if (os.path.exists(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_x.txt')):
        os.remove(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_x.txt')
        os.remove(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_y.txt')

    # 初始x，y
    x = np.zeros(times)
    y = np.zeros(times)

    x[0] = x_init
    y[0] = y_init

    # 初始控制参数
    w = np.arange(0.1, 4, 0.1)
    u = np.arange(0.1, 4, 0.1)

    # 迭代sine_Logistic映射
    for u_value in range(u.size):
        for w_value in range(w.size):
            for n in range(1, times):
                x[n] = np.mod(u[u_value] * np.sin((np.pi * np.pi * y[n - 1]) / (1 - x[n - 1])),1)
                y[n] = np.mod(w[w_value] * np.sin((np.pi * np.pi) / (x[n - 1] * y[n - 1])),1)

                result_x = np.array([u[u_value], w[w_value], x[n]])
                result_y = np.array([u[u_value], w[w_value], y[n]])
                # 存储结果
                with open(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_x.txt', 'a+') as file:
                    np.savetxt(file, result_x, fmt='%f', delimiter=',')
                with open(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_y.txt', 'a+') as file:
                    np.savetxt(file, result_y, fmt='%f', delimiter=',')
    # 装载数据
    x = np.loadtxt(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_x.txt')
    x = x.reshape(int(x.size / 3), 3)

    y = np.loadtxt(r'../experiment/2d_sine_logistic/2d_sine_logistic_biff_y.txt')
    y = y.reshape(int(y.size / 3), 3)
    # 使用scatter画x的三维图
    # fig = plt.figure()
    # ax = Axes3D(fig)
    # fig.add_axes(ax)
    # ax.scatter(x[:, 0], x[:, 1], x[:, 2], )
    # 使用scatter画y的三维图
    fig = plt.figure()
    ax = Axes3D(fig)
    fig.add_axes(ax)
    ax.scatter(y[:, 0], y[:, 1], y[:, 2])

    # 设置坐标系
    ax.set_zlabel('y', fontdict={'size': 15, 'color': 'red'})
    ax.set_ylabel('u', fontdict={'size': 15, 'color': 'red'})
    ax.set_xlabel('w', fontdict={'size': 15, 'color': 'red'})
    ax.set_title("2D Sine-Logistic")

    #显示图片
    plt.show()




def Logistic(x,n):
    for i in range(n):
        y = 4 * x * (1 - x)
        x = y
    return x


def LE_calculate(n = 5000):
    sum_Lambda1 = 0
    a = 0.123456789
    # 使用符号方式求解
    x = symbols("x")
    f_mat = Matrix([4 * x * (1 - x)])
    # 求解雅各比矩阵
    jacobi_mat = f_mat.jacobian([x])  # 带变量的雅各比矩阵形式是固定的
    a = Logistic(a, 5001)  # 先迭代5000次，消除初始影响.以第5001次的值作为初始值
    U1 = Matrix([1])  # 初始列向量
    for i in range(n):
        J = jacobi_mat.subs(x,a)  # 将变量替换为当前迭代值，得到当前的雅各比矩阵（数字）
        column_vector1 = U1#初始列向量为上一次的U1
        vector1 = J * column_vector1  # 初始列向量乘上雅各比矩阵之后得到的向量
        V1 = vector1  # 将vector1复制给V1
        U1 = V1 / (V1.norm(2))  # 向量U1等于向量V1除以向量V1的模(2范数)
        Lambda1 = ln(V1.norm(2))
        sum_Lambda1 = sum_Lambda1 + Lambda1
        a= Logistic(a,1)#进行下一次迭代
    LE1=sum_Lambda1/n
    print(LE1)

def _2D_sine_logistic(u,w,init_x,init_y,times):
    '''
    :param times: 循环次数
    :param u: 系统控制参数 [0,4]
    :param w: 系统控制参数 大于等于2
    :param init_x: 系统初始值
    :param init_y: 系统初始值
    :return: 2Dsine_Logistic迭代结果
    '''

    #赋值初始值
    x = init_x
    y = init_y
    for n in range(1, times):
        x1 = u * np.sin(np.pi*np.pi * y/(1-x))
        y1 = w * np.sin(np.pi*np.pi / (x*y))
        x = x1
        y = y1
    return x, y

def LE_calculate(u , w, n = 5000):
    sum_Lambda1 = 0
    sum_Lambda2 = 0
    # u = 0.1
    # w = 0.1
    a = 0.123456789
    b = 0.123456789
    # 使用符号方式求解
    x, y = symbols("x,y")
    f_mat = Matrix([u * sin(np.pi * np.pi * y/(1 - x)), w * sin(np.pi * np.pi/(x * y))])
    # 求解雅各比矩阵
    jacobi_mat = f_mat.jacobian([x, y])  # 带变量的雅各比矩阵形式是固定的
    a, b = _2D_sine_logistic(u, w, a, b, n+1)  # 先迭代5000次，消除初始影响.以第5001次的值作为初始值
    U1 = Matrix([1, 0])  # 初始列向量
    U2 = Matrix([0, 1])
    for i in range(n):
        J = jacobi_mat.subs({x: a, y: b})  # 将变量替换为当前迭代值，得到当前的雅各比矩阵（数字）
        column_vector1 = U1  # 初始列向量为上一次的U1和U2
        column_vector2 = U2
        vector1 = J * column_vector1  # 初始列向量乘上雅各比矩阵之后得到的向量
        vector2 = J * column_vector2
        V1 = vector1  # 将vector1复制给V1
        U1 = V1 / (V1.norm(2))  # 向量U1等于向量V1除以向量V1的模(2范数)
        V2 = vector2 - (vector2.dot(U1)) * U1  # dot为点乘(內积)
        U2 = V2 / (V2.norm(2))
        Lambda1 = ln(V1.norm(2))
        Lambda2 = ln(V2.norm(2))
        sum_Lambda1 = sum_Lambda1 + Lambda1
        sum_Lambda2 = sum_Lambda2 + Lambda2
        a, b = _2D_sine_logistic(u, w, a, b, 1)  # 进行下一次迭代
    LE1 = sum_Lambda1 / n
    LE2 = sum_Lambda2 / n
    # print('LE of X:{:.5f}'.format(LE1))
    # print('LE of Y:{:.5f}'.format(LE2))
    return LE1,LE2


if __name__ == '__main__':

    u = np.arange(0.1, 4.1, 0.1)
    w = np.arange(0.1, 4.1, 0.1)
    # u,w = np.meshgrid(u,w)
    # print(u)
    length = u.shape
    le1_copy = np.zeros([length[0], length[0]])
    le2_copy = np.zeros([length[0], length[0]])
    print(le1_copy.shape)
    for i in range(length[0]):
        for j in range(length[0]):
            le1,le2 = LE_calculate(u[i], w[j],n=100)
            print('u is :{:.5f}, LE of X:{:.5f}'.format(u[i], le1))
            print('w is :{:.5f}, LE of Y:{:.5f}'.format(w[j], le2))
            le1_copy[i,j] = le1
            le2_copy[i,j] = le2

    np.save('./le1.npy', le1_copy)
    np.save('./le2.npy', le2_copy)


    # # 使用scatter画y的三维图
    # fig = plt.figure()
    # ax = Axes3D(fig)
    # fig.add_axes(ax)
    #
    # u_index = np.zeros([length[0], length[0]])
    # w_index = np.zeros([length[0], length[0]])
    # le1_copy = np.load('le1.npy')
    # le2_copy = np.load('le2.npy')
    # for i in range(length[0]):
    #     for j in range(length[0]):
    #         u_index[i, j] = u[j]
    #         w_index[i, j] = w[j]
    #
    #
    # u_index = u_index.reshape(-1)
    # w_index = w_index.reshape(-1)
    # le1_copy = le1_copy.reshape(-1)
    # le2_copy = le2_copy.reshape(-1)
    #
    # ax.plot_surface(u_index, w_index, le1_copy)
    # # ax.scatter(u_index, w_index, le1_copy, 'r.')
    # # 显示图片
    # plt.show()


